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Introduction - Many characteristics of the artificial satellite 
orbits useful to geodesy are similar to those of the planetary and 
natural satellite orbits which have long been the object of attention 
of celestial mechanics: slowly varying Keplerian ellipses of mooerate 

fP 

eccentricity.. Hence many of the classical methods can be applied. 

> Characteristics peculiar to close satellites* and hence requiring 
^ rSN special attention, are the capability of a wide variety of inclinations, 
^ the dominant perturbation by the earth’s oblateness; the many other 

/small perturbations by variations of the earth’s gravitational field* 

vv v* 

xwhich are the chief reason why satellites are of geodetic interest, the 

x ©importance of atmospheric drag; and the irregular distribution of obser- 

o tc 
« o 

>t xvations around the orbit. These characteristics, together with the 

availability of large high-speed computers, create some interesting new 
problems for satellite geodesy, despite the basic similarity to classical 


Celestial Mechanics - For purely gravitational effects on the orbit 
of a particle, the differential equations can be written: 


'i - 


when £ is the position vector, V is a scalar potential, and V denotes the 
gradient. The three second order equations (i) can easily be converted 
to six first order equations by taking as variables in addition to the 
three position components V the three velocity components r. Then given a 
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law for V and position and velocity at a particular time, V and V , the 
equations of motion can be readily solved by numerical integration. 

However, other developments may be more conducive to insight or more 
economical of computer time. In the first approximation for the earth, 

V ^ kM/r (2) 


where k is the gravitational constant, M is the earth's mass, and r is 
the radius vector. Use of (2) in (l) and conversion of (l) to spherical 
coordinates obtains as a first integral 

r 2 f = h, (3) 

where f is the rate of angular motion about the origin in the plane defined 
by the position and velocity vectors of the satellite, and h is a constant. 
(3) is known as the integral of areas. Replacement of r by its inverse and 
use of (3) to eliminate dt in favor of df results in an equation which can be 
integrated to obtain 1/r in terms of f: 

1_ = 1 + e cos f 

r a( 1-e 2 ) (4) 

the equation for the radius vector of an ellipse with a semi 'major axis a, 
an eccentricity e, the origin at a focus, and f, called the true anomaly, 
measured from the point of closest approach to the focus. We also get 

h = [kHa (l-e 2 )J 


( 5 ) 
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To completely specify the position of the particle, in addition to a,e, 
and f three more variables are required to express the orientation cf the 
ellipse in the plane of motion: <0, the argument of perigee; and to express 
the orientation of the plane with respect to an inertial reference: i, the 

inclination, and il, the longitude of the node. See Figure 1 for the definition 
of these variables., vkv\V t<k*Vi««J[ W\lV\ "Hi- 

To integrate (3) to obtain the position in terms of time, we replace 
the true anomaly f by the eccentric anomaly E as follows: define the major 

axis as the abscissa; take a circle tangent to the ellipse at the ends of the 
major axis; and define the eccentric anomaly E as the angle subtended at the 
center of the circle by the point on the circle whose abscissa is the same as 
the point on the ellipse of true anomaly f. Then 

T» ■= a(l - e cos E) (6) 


Differentiating (4) and (6) and equating the two expressions for dr, df 

can be eliminated for dE and a form of (3) is obtained which can be integrated 

to: 


E 


e sin E 


^ .' 3/2 



= 71 (t - t ) = M 

0 


CO 




O : 


where t o is the time of passing the point of closest approach 
called perigee'^C?) is called Kepler's equation. 

0f\, Vs uiU I \W J (Ht\i M i*> tdhJi 4 \Sl MM 
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For cases where the potential V in (l) does not have the central form 
(2), but is instead: 


V 



R, 


(8) 


where R is known as the disturbing function, the Keplerian ellipse and its 
orientation can still be regarded as a coordinate system: at any instant the 

situation of a satellite in inertially fixed geocentric coordinates can be 
expressed by the six numbers |a,e,i ,M,y),Q^ instead of the six numbers 

the Keplerian elements, the six differential equations in terms of^ir.i^can be 
converted, after considerable algebra, into: 


y,z,x 


,y,zf. By using the partial derivatives of r, r with respect 

Jr Wi 


to 


da _ 2_ 

dt ” na cfi 

de _ 1_ - el £R - (l - e 2 A » 

dt ~ na*e 2M na 2 e dm 

i 

dm cos i i „ fiR + (l - e 2 ) . £R 

dt " na*(l-e)* sin i di na 2 e de 

di _ cos i , dR - 1 . . SR 

dt ” na^l-e 2 )? sin i du na 2 (l-e 2 )8 sin i Sfi 


veo 


dCl . 1 i . SR 

dt ” nd 2 ( l-e a )*sin i di 

dM n _ 1 - e 2 . » . 2 . SR 

dt = " Na^e de " na da 

J 

This form of the equations of motion is known as the Lagrangian planetary 
equations. 


The Earth's Gravitational Field . To use equation ( 9 ), the potential 
V, including the disturbing function R, must be expressed in terms of Keplerian 
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elements. Because of tne strong smoothing effects of the extrapolation to 
altitude and the integration to obtain satellite position, the most convenient 
representation of the Earth's gravitational field for application to orbits is 
in spherical harmonics (see "Use of Spherical Harmonics in Physical Geodesy" 
by U. A. Uotila and "The Figure of the Earth" by W. M. Kaula) : 

WM ® n n 

V = — (.1 + E £ ( — ) P (sin Y) (C cos m\ + S sin m\) ] (lO) 

r n=2 ra=0 r nm nm nm 

where a« is the equatorial radius, Y is the geocentic latitude, P (sin Y) 

^ nm 

is the Legendre Associated Polynomial, \ is the longitude, and „ S^ m are 

independent coefficients. The interest of geodesy in close satellite orbits 

is to use their perturbations to determine these coefficients C , S By a 
r nm nm 

purely geometrical transformation, treating the Keplerian ellipse as an alter- 
native to spherical polar coordinates, there is obtained from ( 10) a disturbing 
function (Kaula, 1%1) : 

* = f y & y i t (i> i 6 (.), 

& a /— » L-i nmp Z_ j np cj 

n=2 m=0 P— 0 q--® 

- . -n-m even 

X Li-s" 1 } COS {( n *2p) U) + (n-2p+q)H + m(Q-6)} 
nm n ~m odd 

S nm . n-m even 

+ |c nm | sin { ( n - 2p ) uc + (n-2p+q) H + ra(Q-0) } j 

n-m odd * * 


where 6 is the Greenwich Sidereal Time, and 
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, r (2n-2t) j 

^nmp'*' = L, t I (n-t) I (n-m- 2 t) l 2 ^ n ” 2 t 


. n-m-2t . 
sin i x 


t 

m 

X 

s=0 


iO co®* 1 i( n '”c 2t+s )(p"t!c) ( ' 1)C 


-k 


( 12 ) 


In (12), k is the integer part of (n-m)/2; t is summed from 0 to the 
lesser of p or k, and c is summed over all values making the binomial coef- 
ficients non-zero* G^^e) is always 0(e^ )• The general expression for 
G^^'e) is rather complicated* There is, however, a simpler form for the 
more important long period terms for which M is absent from the argument; 
i.e*, for which q (2p-n): 

, \ 1 \ f n-l ./2d+n-2p , ^ev 2d+n *' 2pl 

G np(2p-n) 6 = ( l-e 2 ) n ”^ L \2d+n-2pVv d }V 

d=0 


(13) 


where p* = p for p£ n/2, and 

p* = p for p^ n/2. 

The form of ( 1 1 ) indicates that the elements a.e.i appear in every term 
of the disturbing function R, but that by an appropriate combination of sub* 
T’-i nts n,m,p,q the elements u.',M,0 and the sidereal time 9 may be absent from 
certain terms. One such combination is nmp^ = 2010, for which the coefficient 
C 20 is the dominant coefficient in R, expressing the earth’s oblateness. Use 
of this term as the disturbing function in the equations of motion ( 0 ) results 
in no perturbation of a,e, or i, but in constant changes in uj, 0, and M: 
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(f*)' 2 (1-5 cos 2 i) « 3.35 (5 cos 2 i - L) degrees/day 

= ^ 77 ^ 2 x 2 (^ e -) 2 cos i rs - 6.70 cos i degrees/day 
at i *e ; a 

= n - ^ 2)rj A) 2 (3 cos 2 i - l) « 14 . 3 '' + .0033 (3 cos 2 i - 1) revs/day 


where (12) and ( 13) have been used to obtain F<oi(i) and G 2 io(e), and the 
numerical examples are based on a C 20 of -.000434 and typical orbit specifications 
of .01 eccentricity and 1.12 a e semi -major axis.. 

Rates of secular motion of the perigee and node as specified by (14) are 
the dominant characteristic observed in satellites which are high enough not to 
be in a state of rar ui decay due to atmospheric drag. This characteristic is so 
dominant that for the effect of other terms in the earth's field, as given by 
(ll), the integration can be performed assuming that a,e» and i are con cra nr and 
u>>vl, and M have a constant rate of change on the right hand side of (?>. For 
example, the perturbation of the node by a terra of subscripts n^p,© 1 



c n m even 

I { ” m } sin {(n-2p)uj + (n-2pHq) M + m(Q-0)} 
L nm n-m odd 


r Snm -) 
■ICnm 3 


n-m even 
cos 

n • m odd 


((n-2p) u) + 


(n-2p-tq)M 



(15) 
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A particular term which has been the object of considerable orbit analysis 
is A£! 22 io, the semi -daily perturbation of the node by the equatorial ellipti- 
city: 


AQ221 0 


k Ma|(3F 2 ail3i)G2i o 
2na^( l-e 2 )®sin i(fi-§) 


[C 22 sin 2(Cl-0) 


S 22 cos 2( Cl- 0 ) ] 


= [° 22 sin 2(Q * 9) ' S22 cos 2(CJ - 0 ^ (l6 ) 

• • • 

The appearance of the rate {(n-2p)<L + (n-2p4q)M •+• ra(Q-9)) in the denomina- 
tor of (15) and the rates given by (14) indicate that the largest perturbations 
should be those for which (n-2p+q) and m both vanish: i.e., long period 01 
secular effects of the zonal harmonics C no , and that the largest effect of a 
tesseral harmonic C nm , S nro (for which m^O) will normally be one for which 
(n-2p+q) vanishes. 

Non-linear Perturbations. Th description of the orbit as a Keplerian 
ellipse described by equations (6) and (7) plus secular motions C 14) and linear 
periodic perturbations such as (15) is inadequate for geodetic purposes. Because 
the oblateness C 2 o is about 1000 times larger than the other coefficients C nm , 
S nm , there are to be expected non-linear effects of coefficient C 2 o which are of 
comparable magnitude to the linear effects of coefficients C nm , S nra . As in any 
non-linear problem, there are many ways to proceed. Some theories develop (9) 
or similar equations to the next higher order by substituting the linear approxi- 
mation on the right (e.g., Kozai, 1759b). The algebra entailed in such a develop 

t 

ment is considerable, however. An alternative frequently adopted is to transform 
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(9) to a simpler form by a change of variables; the most popular such set are 


the Delaunnay variables'. 

L = (na)*, /, = M 

G - [^a(l-e 2 )]^, g = (D 

H = [>a( l -e 2 ) J^Coslj h = Q 

If we set 


(17) 


(18) 


known as the force funccion, then (9) becomes: 

dL _ 3F dt _ SF 

dt = bV dt ~ " $L 


dG J dj JF 
dt = &U* dt = ' SG 


dH _ jJF dh = . £F 

dt ~ SCI* dt 2H 


(19) 


Normally a closed solution to (19) with F as in (13) cannot be found. The 
solution customarily adopted is to take the known solution for another force 
function F^lose to F and express the solution for F as Tayior series for the 
small differences L-L* , G-G* , etc. Continuity conditions require that the terms 
of these Taylor series all be expressible as derivatives of a single scalar 
function, called the determining function. The solution for F* is called the 
intermediate orbit. Different theories usually differ in their selection of the 
intermediate orbit. A theory which has a simple intermediary will have compli- 
cations in developing the determining function and its derivatives (e.g., Brouwer, 
1959; Kozai , 1962); while a theory which has an intermediary which accounts for 
nearly all the motion will entail considerable algebra in expressing position and 
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velocity in the intermediary (e.g. , Vinti, 1959. 1961; Izsak, 1960). 

Some theories have intermediaries which are not defined by a force 
function F* , but rather by geometrical conditions, such as those of Merson 
(i960 and Musen (1959, I960. The theory of Musen (1959) also differs in that 
it is designed for solution by numerical iteration, rather than by analytical 
series development. 

Most standard close satellite theories break down for orbits which have 
an eccentricity or inclination around zero. The remedy is to use orbital 
elements which do not require a distinct perigee or node for their precise 
numerical definition. 

Another cause of breakdown is that for some combination of subscripts 
n,m,p,q the rate appearing in the denominator of the perturbation, such as 
(15)> approaches zero: 

(n-2p) ub + (n-2p+q)M + m(6-0) « 0, (20) 

so that a resonant situation occurs. The one which has received the most atten- 
tion is that of the critical inclination, for which & is zero (e.g., Izsak, 1%3). 
Resonant orbits of more interest to geodesy are those of synchronous communication 
satellites, for which (<L + ft + Cl - 6) is zero. These orbits will resonate with 
any term in the gravitational field for which (n-m) is even, since it will have 
the near-zero rate for a subscript p of (n-m)/2 and q zero. However, the semi- 
major axis of such an orbit is so large- -about 6.6 earth radii-- that it probably 
would be a sensitive determiner of only the lowest degree such harmonic: C 22 , S 22 . 
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The stable point for the resonance will be the longitude of the minimum for the 
equatorial eillpticity, [tan- 1 (S 22 /C 22 ) ± tt]/2. In the first approximation, 
the solution is similar to that for a pendulum, with the longitude of the 
satellite expressed by an elliptic integral. If the initial longitude is 
sufficiently close to the stability point, the longitude will librate about the 
point; if it is outside these limits, the satellite will drift. The rate of 
motion will be proportionate to (C?i + S 2 2 )^. See Musen and Bailie (1962) and 
Morando (1963)- 

Other Perturbations. The most important other effect on close satellites 
of geodetic interest is atmospheric drag. Because of the much greater density 
in vicinity of perigee, the principal effect of drag is a slow decrease in the 
eccentricity and the semi -major axis. The rate of this decrease will vary with 
the rotation of the perigee with respect to the bulge of the atmosphere toward 
the sun, and with the intensity of the solar flux. In addition, there will be 
an effect on the rates of motion of the node and perigee due to interaction 
between the drag and the earth's oblateness, which must be considered in ana- 
lyzing secular changes to determine the even degree zonal harmonics. All these 
variations make it impracticable to develop an analytical theory of drag. 
Numerical integration or numerical harmonic analysis is necessary to apply any 
realistic atmospheric model. Furthermore, there are incalculable irregular 
variations with respect to even the best of atmospheric models; as a consequence 
arbitrary polynomials are more frequently used. For geodetic satellites it is 
desirable to keep the perigee high enough and area-to-mass ratio low enough, so 
that irregular variations whose wave lengths are not long with respect ;o the 
spacing of observations are negligible. For perigees above 300 or 1000 kilo- 
meter altitude and area-to-mass ratios less than 0.2 cm 2 gm- 1 , the amplitude of 
irregular variations of less than monthly period are generally less than 5 x 10- 5 
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